%% exactly result
clear all
angle = [0.1:0.1:atan(2) - 0.05];
liquid_volume = 100;
s = sin(angle);
c = cos(angle);
t = s./ c;
r = 5;
h = 20;
ex_level = 2 * r * s + (liquid_volume * c - 2 * r * r * s) / (2 * r);
%% monte carlo result
x = [0 20];
y = 1 + [4 4];
[bottle_mass_center bottle_shape] = bottle_model(x, y);
angle_size = size(angle);
angle_end = angle_size(2);
mc_level = [];
for i = 1:angle_end
    rotated_bottle_shape = rotation_matrix(-angle(i)) * bottle_shape;
    mc_level(i) = liquid_level(liquid_volume, rotated_bottle_shape);
   
    plot(rotated_bottle_shape(1,:), rotated_bottle_shape(2,:), 'k-');
    hold on
    plot([min(rotated_bottle_shape(1,:)) max(rotated_bottle_shape(1,:))],...
         [mc_level(i) mc_level(i)], 'b-');
    plot([min(rotated_bottle_shape(1,:)) max(rotated_bottle_shape(1,:))],...
         [ex_level(i) ex_level(i)], 'r-'); 
    ylim([0 20]);
    xlim([-15 20]); 
    axis equal
    M(i) = getframe;
    hold off
end
movie(M, 100);
figure;
plot(angle, ex_level, 'rx-');
hold on
plot(angle, mc_level, 'bo-');
legend(3, 'Exactly', 'Monte Carlo');
ex_area = 2 * r * r * t + (ex_level - 2 * r * s) * 2 * r./ c;
mc_area = 2 * r * r * t + (mc_level - 2 * r * s) * 2 * r./ c;
figure;
plot(angle, ex_area, 'r');
hold on
plot(angle, mc_area, 'b');
legend(3, 'Exactly', 'Monte Carlo');